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1.  INTRODUCTION 

Networks  of  queues  occur  frequently  in  diverse  applications.  In  particular,  they  are  widely  used 
in  studies  of  computer  and  communication  system  performance  as  models  for  the  interactions  among 
system  resources.  This  paper  deals  with  mathematical  and  statistical  methods  for  discrete  event 
simulation  of  networks  of  queues.  The  emphasis  is  on  methods  for  the  estimation  of  general  characteris¬ 
tics  of  "passage  times"  in  closed  networks.  Informally,  a  passage  time  is  the  time  for  a  job  to  traverse  a 
portion  of  a  network.  Such  quantities,  calculated  as  random  sums  of  queueing  times,  are  important  in 
computer  and  communication  system  models  where  they  represent  job  response  times. 

The  most  obvious  methodological  advantage  of  simulation  is  that  in  principle  it  is  applicable  to 
stochastic  systems  of  arbitrary  complexity.  In  practice,  however,  it  is  often  a  decidedly  nontrivial  matter 
to  obtain  from  a  simulation  information  which  is  both  useful  and  accurate,  and  to  obtain  it  in  an 
efficient  manner.  These  difficulties  arise  primarily  from  the  inherent  variability  in  a  stochastic  system, 
and  it  is  necessary  to  seek  theoretically  sound  and  computationally  efficient  methods  for  carrying  out  the 
simulation.  It  is  fundamental  for  simulation,  since  results  are  based  on  observation  of  a  stochastic 
system,  that  some  assessment  of  the  precision  of  results  be  provided.  Assessing  the  statistical  precision 
of  a  point  estimate  requires  careful  design  of  the  simulation  experiments  and  analysis  of  the  simulation 
output.  In  general,  the  desired  statistical  precision  takes  the  form  of  a  confidence  interval. 

The  estimation  methods  described  here  for  passage  times  in  networks  of  queues  use  the 
regenerative  method  [2],  [3],  [5]  for  analysis  of  simulation  output.  Based  on  a  single  simulation  run. 
they  provide  (strongly  consistent)  point  estimates  and  (asymptotically)  valid  confidence  intervals  for 
general  characteristics  of  "steady  state"  (limiting)  passage  times. 

The  notion  of  a  distinguished  "marked  job"  is  fundamental  to  the  method  for  estimation  of 
passage  times  described  in  Section  3.  We  arbitrarily  select  a  job  to  serve  as  the  marked  job.  and 
measure  its  passage  times  in  the  course  of  the  simulation.  In  developing  the  marked  job  method,  we 
take  as  the  state  vector  of  the  network  of  queues  the  usual  numbers-in-queue,  stages-of-service  state 
vector  augmented  by  information  sufficient  to  track  the  marked  job.  A  key  step  in  the  development  of 
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the  method  is  the  definition  of  a  related  continuous  time  stochastic  process.  The  starts  of  passage  times 
for  the  marked  job  are  the  successive  entrances  of  this  process  to  a  fixed  set  of  states,  and  the  termina¬ 
tions  of  such  passage  times  are  the  successive  entrances  to  another  fixed  set  of  states.  We  develop  in 
this  setting  a  ratio  formula  from  which  point  estimates  and  confidence  intervals  can  be  obtained  for 
quantities  associated  with  the  limiting  passage  time. 

In  Section  4  we  give  a  second  estimation  method  for  passage  times.  With  this  decomposition 
method,  observed  passage  times  for  all  of  the  jobs  enter  into  the  construction  of  the  point  and  interval 
estimates.  This  estimation  method  applies  to  the  restricted  but  practically  important  class  of  passage 
times  which  are  not  complete  circuits  in  a  network. 

2.  NETWORKS  OF  QUEUES  AND  PASSAGE  TIMES 

We  deal  with  closed  networks  of  queues  having  a  finite  number  of  jobs  (customers),  N,  a  finite 
number  of  service  centers,  s,  and  a  finite  number  of  (mutually  exclusive)  job  classes,  c.  At  every  epoch 
of  continuous  time  each  job  is  in  exactly  one  job  class,  but  jobs  may  change  class  as  they  traverse  the 
network.  Upon  completion  of  service  at  center  i  a  job  of  class  j  goes  to  center  k  and  changes  to  class  1 
with  probability  pjj  k|,  where  is  a  given  irreducible  Markov  matrix.  At  each  service  center 

jobs  queue  and  receive  service  according  to  a  fixed  priority  scheme  among  classes;  the  priority  scheme 
may  differ  from  center  to  center.  Within  a  class  at  a  center,  jobs  receive  service  according  to  a  fixed 
queue  service  discipline,  e.g.,  first-come,  first-served  (FCFS).  Note  that  in  accordance  with  the  matrix 
P.  some  centers  may  never  see  jobs  of  certain  classes.  Only  one  job  can  receive  service  at  a  center  at  a 
time.  i.e..  the  centers  are  single  servers.  According  to  a  fixed  procedure  for  each  center,  a  job  in  service 
may  or  may  not  be  preempted  if  another  job  of  higher  priority  joins  the  queue  at  the  center. 

Probabilistic  assumptions 

The  marked  job  method  discussed  in  Section  3  for  passage  time  simulation  applies  to  networks 
of  this  kind  in  which  all  service  times  in  the  network  are  mutually  independent,  and  at  a  center  have  a 
distribution  with  a  Cox-phase  representation  [1],  i.e..  consisting  of  a  sequence  of  exponential  stages;  see 
Figure  1.  We  permit  parameters  of  the  service  time  distribution  to  depend  on  the  service  center,  the 
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class  of  job  being  served,  and  the  "state"  of  the  entire  network. 

Each  service  time  distribution  has  its  own  finite  number  of  stages,  say  n.  Realization  of  a 
service  time  is  as  a  sum  of  a  random  number  (£l)  of  exponentially  distributed  times.  Within  the  jth 
stage  (l<jsn-l).  an  amount  of  service  exponentially  distributed  with  parameter  Xj  accrues;  with 
probability  1-bj  service  is  complete,  and  with  probability  bj  additional  service  exponentially  distributed 
with  parameter  XjH  accrues.  The  density  function  of  the  resulting  service  time  has  rational  Laplace 
transform 

f (s)  -  2  (l-b.)a  n  Xk/(Xk-s)  . 

j-i  k-l 

where  a,-l  and  for  l<jsn,  aj—bj ...bj.i .  The  class  of  density  functions  having  Laplace  transforms  of  this 
form  includes  hyperexponential  and  mixtures  of  Erlang  densities;  the  Appendix  of  [4],  Note  that  we 
exclude  the  case  of  zero  service  times  occurring  with  positive  probability. 

State  vector  definition 

In  order  to  characterize  the  state  of  the  network  at  time  t,  we  let  Sj(t)  denote  the  class  of  the  job 

receiving  service  at  center  i  at  time  t,  where  i— 1.2 . s;  by  convention  Sf(t)-0  if  at  time  t  there  are  no 

jobs  at  center  i.  The  classes  of  jobs  serviced  at  center  i  ordered  by  decreasing  priority  are  j,(i).  j;(i) . 

jkji)(i).  elements  of  the  set  { 1,2 . c}.  Let  Cj)<i>(t) . Cj  (i)(t)  denote  the  number  of  jobs  in  queue  at 

time  t  of  the  various  classes  of  jobs  serviced  at  center  i.  i-1,2 . s.  For  the  queue  length  of  jobs  of 

various  classes  at  the  several  centers,  these  state  variables  (together  with  the  stages-of-service)  would 
suffice.  They  are  not  adequate,  however,  to  deal  with  general  characteristics  of  passage  times.  An 
apparently  minimal  state  vector  augmentation  is  based  on  the  concept  of  a  marked  job.  The  idea  is  to 
keep  track  of  the  position  of  an  arbitrarily  chosen  marked  job.  and  to  measure  its  passage  times  during 
the  simulation.  It  is  convenient  to  think  of  the  N  jobs  being  completely  ordered  in  a  linear  stack 


according  to  the  following  scheme.  For  t^O.  we  define  the  vector  Z(t)  by 

z(l>  "  (Cj«i,(,,(t) . CjI<,,^t)’S,(t):-:CJk(,l(5)(t) . ■Cj|‘,»(t).Sl(t)). 


(1) 
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The  linear  job  stack  then  corresponds  to  the  order  of  components  in  the  vector  Z(t)  after  ignoring  any 
zero  components.  Within  a  class  at  a  particular  service  center,  jobs  waiting  appear  in  the  linear  stack  in 
FCFS  order,  i.e..  in  order  of  their  arrival  at  the  center,  the  latest  to  arrive  being  closest  to  the  top  of  the 
stack.  We  denote  by  N(t)  the  position  (from  the  top)  of  the  marked  job  in  this  linear  stack. 

W  ith  each  job  in  the  network,  we  associate  a  stage  of  service  as  follows.  A  job  in  service  is  in  a 
particular  stage  of  its  service  time  distribution  at  that  center;  for  such  a  job.  this  is  the  associated  stage. 
For  a  job  in  queue  at  a  center,  the  associated  stage  is  that  stage  of  service  which  is  to  be  rendered  when 
the  job  next  receives  service;  typically  this  is  the  first  stage  of  service,  but  may  be  a  subsequent  stage  if 
the  job  has  been  preempted.  For  t^O,  we  define  the  vector  U(t)  by 

U(t)  -  (U,(t) . UN(t))  , 

where  Uj(t)  is  the  stage  of  service  associated  with  the  jth  job  in  the  linear  stack  of  jobs,  and  take  as  the 
state  vector  of  the  network  of  queues  the  vector 

X(t)  -  (Z(t).N(t),U(t))  .  (2) 

For  any  service  center  i  that  sees  only  one  class  of  job  (k(i)-l),  it  is  possible  to  simplify  the  state  vector 
by  replacing  Ck  (0(t).  Sj(t)  by  Qj(t),  the  total  number  of  jobs  at  center  i. 

Definition  of  passage  times 

Given  a  particular  closed  network  of  queues,  we  must  specify  the  passage  time  of  interest.  This 
can  be  done  in  terms  of  the  the  marked  job  by  means  of  four  subsets  (A,.  A:.  B,.  B-0  of  the  state  space. 
E.  of  the  stochastic  process  X- { X(t)  U:0 ) .  The  sets  A,.  A,  [respectively  B,.  B,]  jointly  define  the 
random  times  at  which  passage  times  for  the  marked  job  start  [respectively  terminate].  The  sets  A,.  A:, 
B,.  B;  in  effect  determine  when  to  start  and  stop  the  clock  measuring  a  particular  passage  time  of  the 
marked  job. 

It  is  convenient  to  introduce  the  jump  times  [rn:n20J  of  the  process  X.  For  k.nal.  we  require 
that  the  sets  A,.  A2,  B,  and  B;  satisfy; 

if  X(rn.,)«A,.  X(rn)«  A:.  X(Tn.,.k)«  A,  and  Xfr,,^)*  A,, 
then  X(rn.|.m)«B|  and  X(rn_m)*B:  for  some  0<msk; 
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and 

if  X(rn.,)tB|.  X(rn)tB:.  X(Tn.,_k)<B,  and  X(Tn_k)eB2. 
then  X(rn.,,m)«  A;  and  X(rn.m)«A:  for  some  Osm<k; 

These  conditions  ensure  that  the  start  and  termination  times  for  the  specified  passage  time  strictly 
alternate. 

In  terms  of  the  sets  A|,  A2.  Bj,  and  B:.  we  define  two  sequences  of  random  times.  and 

{Tj.jsl  J.  where  Sj.t  is  the  start  time  of  the  jth  passage  time  for  the  marked  job  and  Tj  is  the  termina¬ 
tion  time  of  this  jth  passage  time.  Assuming  that  the  initial  state  of  the  process  X  is  such  that  a  passage 
time  for  the  marked  job  begins  at  t-0,  let 

S0  "  0 

Sj  "  inf { rn£Tj:X(rn) t  Aj.  X(Tn.,)t  A,  J  ,  jal 
and 

Tj  -  inf{Tn>Sj.,;X(rn)tB2<  X(Tn.,)tB| }  ,  j^l  . 

Then  the  jth  passage  time  for  the  marked  job  is  simply  Pj-Tj-Sj.  ],  jal.  For  passage  times  that  arc 
complete  circuits  in  the  network,  A^B,  and  A2-B2;  consequently  Sj— Tj  for  all  jsl. 

3.  THE  MARKED  JOB  METHOD 

Estimation  of  general  characteristics  of  passage  times  by  simulation  of  a  closed  network  of 
queues  can  be  based  on  the  measurement  of  passage  times  for  a  typical  job.  the  marked  job  discussed 
above.  It  is  intuitively  clear  (and  is  shown  in  the  Appendix  of  [6])  that  the  sequence  of  passage  times 
for  any  other  job  (as  well  as  the  sequence  of  passage  times  irrespective  of  job  identity,  in  order  of  start 
or  termination)  converges  in  distribution  to  the  same  random  variable  as  the  sequence  of  passage  times 
for  the  marked  job.  It  follows  that  we  can  estimate  general  characteristics  of  passage  times  in  the 
network  by  simulation  of  the  process  X.  We  use  the  regenerative  method  (applied  to  a  stochastic  process 
defined  in  terms  of  X)  to  obtain  from  a  single  simulation  run  point  and  interval  estimates  for  passage 


time  characteristics. 
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We  denote  by  Xn,  naO.  the  state  of  the  system  when  the  (n-1  tst  passage  time  of  the  marked  job 
starts.  For  ja|.  let  Pj  be  the  jth  passage  time  for  the  marked  job  and  take  the  quantity  of  interest  in  the 
simulation  to  be 

r(f)  “  E  j  f(P)  j  ,  (3) 

where  f  is  a  real-valued  function.  The  quantity  P  is  the  limiting  passage  time.  i.e.. 

lim  P !  Pn=tx  1  -  P{P*x|  . 

n-oe 

for  all  points  x  at  which  the  right  hand  side  of  equation  is  continuous.  For  example,  to  estimate  EjP|. 
we  take  f(pi-p;  to  estimate  PfPstj.  we  take  f(p)-l [0 t](p).  where  1  [0.tj  is  the  indicator  function  of  the 
set  [O.t] . 

Algorithm  Marked  job  method 

1.  To  serve  as  a  return  state,  select  a  state  of  the  system.  i0,  at  which  a  passage  time  of  the  marked  job 
starts.  Begin  the  simulation  with  X(0)-i0. 

2.  Carry  out  the  simulation  of  X  for  a  fixed  number,  n,  of  cycles  (having  random  length)  defined  by  the 
successive  returns  to  the  state  i0. 

3.  In  each  cycle  measure  all  the  passage  times  of  the  marked  job  and  record  these  along  with  the 
number  of  passage  times  for  the  marked  job  in  the  cycle. 

4.  For  kal.  denote  the  number  of  passage  times  observed  in  the  kth  cycle  by  Mk  and  compute  Yk(f), 
the  sum  of  the  quantities  f(Pj)  for  the  passage  times  Pj  in  the  kth  cycle. 

5.  Take  as  a  point  estimate  (based  on  n  cycles)  of  r(f)  the  quantity 

rn(f)  -  Yn(f)/Mn  . 

Yn(f)  -  n*1  1  Yk(f) 

k-l 


where 
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and 

2  Mk. 

k-l 

6.  Take  as  a  100  (l-Z-y)^  confidence  interval  (based  on  n  cycles)  for  r(f)  the  interval 

ln(f)  -  Irntn-Z^^Sn/fMnn'/2).  rn(f)-Z,.,Sn/(Mnn'  :)]  . 

Here  z^-*'1)  1-y),  where  ^t*)  is  the  distribution  function  of  a  standardized  (mean  zero,  variance  one) 
normal  random  variable,  and  sn  is  the  quantity 

sn  “  [S|r2rn(f)si2+(rn):s::]1/* 

where  s, s;;.  and  s12  are  the  usual  unbiased  estimates  for  Var  { Yk(f) } .  Var  {  Mk } ,  and 
Cov{Yk(f).  MkJ.  respectively. 

The  underlying  stochastic  structure 

The  force  of  the  assumptions  of  Cox-phase  service  time  distributions  and  Markovian  routing  in 
the  closed  networks  of  queues  is  that  the  process  X  is  an  irreducible,  positive  recurrent  continuous  time 
Markov  chain  with  a  finite  state  space.  E.  We  let  Xn  denote  the  state  of  the  Markov  chain  X  when  the 
(n-l)st  passage  tim'  of  the  marked  job  begins:  Xn-X(Sn),  naO.  The  process  {Xn:ns0}  is  a  discrete 
time  Markov  chain  (with  state  space  A;)  which  we  assume  to  be  irreducible  and  aperiodic.  Next  we 
observe  that  the  process 

{(X„.Pn-,):n*0} 

is  a  regenerative  process  in  discrete  time.  The  regenerative  property  guarantees  ([11])  that  as  n-co. 

(Xn.Pn.,H(X.P). 

where  —  denotes  weak  convergence  (convergence  in  distribution).  The  random  variables  X  and  P  are. 
respectively,  the  limiting  state  vector  for  the  network  and  the  limiting  passage  time  for  the  marked  job. 

The  goal  of  the  simulation  is  estimation  of  r(f)  -  E[f(P)}.  where  f  is  a  real-valued  measurable  function 
with  domain  R_-[0,f<x>). 

Li  . « 
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We  denote  b>  L(t)  the  last  state  visited  by  the  Markov  chain  X  before  jumping  to  X(t).  and  for 
t^O  define 

V(t)  -  (L(t).X(t))  .  (4) 

The  process  V-jV(t):taO{  is  the  fundamental  stochastic  process  of  the  passage  time  simulation.  This 
process  \  has  a  state  space  F  consisting  of  all  pairs  of  states  (i.j).  i,j<  E.  for  which  a  transition  in  X  from 
state  i  to  state  j  can  occur  with  positive  probability.  Since  X  is  an  irreducible,  positive  recurrent  Markov 
chain,  so  is  V.  We  define  subsets  S  and  T  of  F  according  to 

S-  j  (k,m)  t  F:k  t  A,,  m t  A: ) 

and 

T-  { (k.m)<  F:kt  B,.  m  t  B-. }  (5) 

and  observe  that  the  entrances  of  V  to  S  [resp.  T]  correspond  to  the  starts  [resp.  terminations]  of 
passage  times  for  the  marked  job. 

The  next  step  is  to  select  a  fixed  element  of  S.  which  for  convenience  we  designate  state  0.  We 
set  V(0)-0  and  let  (Vn:naOJ  denote  the  embedded  jump  chain  associated  with  the  continuous  time 
process  V.  We  denote  by  jyn:nal}  the  lengths  in  discrete  time  units  of  the  successive  0-cycles 
(successive  returns  to  the  fixed  state  0)  for  {Vn:ns0j,  and  define  50-0  and  om— y, >m.  msl.  Then 
the  number  of  passage  times  for  the  marked  job  in  the  first  0-cvcle  of  V  is 


and  the  sum  of  the  values  of  the  function  f  for  the  passage  times  for  the  marked  job  in  this  cycle  is 

M, 

Y,  -  I  f(P,)  . 

j- 1 

We  denote  the  analogous  quantities  in  the  kth  0-cycle  by  Mk  and  Yk(ff.  Since  V  is  an  irreducible, 
positive  recurrent  Markov  chain,  it  is  a  regenerative  process,  and  the  pairs  of  random  variables 

{ (Yk(f).Mk):kal }  (6) 


arc  independent  and  identically  distributed  (i.i.d).  Then,  provided  that  P  {  Pt  D(f)  j -0,  where  D(f)  is  the 
set  of  discontinuities  of  the  function  f  in  the  definition  of  r(f).  and  E{  |  f(P)  |  )<oo.  it  follows  that 

r(0  “  E{  Y,(f)  J /E {  M, )  .  (7) 

Given  that  the  random  variables  in  (6)  arc  i.i.d.  and  the  ratio  formula  of  (7),  the  regenerative  method 
applies  and  (from  a  fixed  number  n  of  cycles)  provides  the  strongly  consistent  point  estimate 

rn(f)  -  Vn(f)/Mn 

for  r(f).  The  associated  confidence  interval  is  based  on  the  central  limit  theorem 

n1'  -  j  rn(f)-r(f)  }/(c(f)/E{  M,  ] )— N(0,1 )  .  (8) 

where  <r-(f)  is  the  variance  of  Y , ( f)-r( f)M j  and  N(0.1)  is  a  standardized  (mean  0,  variance  1)  normal 
random  variable.  For  an  alternative  derivation  of  the  marked  job  method  based  on  Markov  renewal 
processes,  see  [7], 

An  Example 

There  are  two  service  centers  in  the  network  of  queues  shown  in  Figure  2.  Upon  completion  of 
exponentially  distributed  a  service  at  center  1,  in  accordance  with  a  binary  random  variable  \f>,  a  job 
joins  the  tail  of  the  queue  in  center  1  (when  i^-l)  or  joins  the  tail  of  the  queue  in  center  2  (when 
After  completion  of  exponentially  distributed  0  service  at  center  2.  the  job  joins  the  tail  of  the  queue  in 
center  1.  Neither  center  1  nor  center  2  service  is  subject  to  interruption.  Assume  that  both  queues  are 
served  according  to  a  first-come,  first-served  (FCFS)  discipline,  and  consider  the  limiting  passage  time  P 
defined  as  the  time  measured  from  entrance  into  the  center  1  queue  upon  completion  of  a  center  2 
service  until  the  job  next  enters  the  center  2  queue. 

In  this  network  there  are  two  classes  of  jobs:  class  1  jobs  at  center  1  and  class  2  jobs  at  center  2. 
Since  each  center  sees  only  only  one  job  class,  by  taking  into  account  the  fixed  number  of  jobs  in  the 
network,  we  can  define  Z(t)  to  be  the  number  of  jobs  waiting  or  in  service  at  center  1  at  time  t.  Then 
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the  process  X-  I  (Z(t).N(t)):tsO  j ,  where  N(t)  is  the  position  of  the  marked  job  in  the  job  stack  at  time 
t.  has  state  space 


E  -  { (i,j):0sisN,  IsjsN}  . 

For  the  passage  time  P.  the  sets  A]  and  A;  defining  the  starts  of  passage  times  for  the  marked  job  are 

A,  -  { (i,N):0si<N  } 
and 

A;  -  { ( i,  1  ):0<isN  J  . 

Similarly.  the  sets  Bj  and  B;  defining  the  terminations  of  the  passage  time  P  are 

B,  -  { (i.i):0<isN  J 
and 

B;  -  { (i*l,i):0<isN )  . 

The  process  V-  j  (L(t).X(t)):taO ) ,  where  L(t)  is  the  last  state  visited  by  the  Markov  chain  X  before 
jumping  to  X(t).  has  state  space 

F-  { (ij.i+l.j-H  ):0^i<N,  lsj<Nj  U  { (i,N,i-l,l):0si<N  } 

U  {(i,j,i-l,j):0<isN,  lsj^N}  U  { (i,i.i,l):l<i^N }  . 

The  subsets  of  F  defining  the  starts  and  terminations  of  passage  times  for  the  marked  job  are 

S  -  { (i.N.i-l,l):0si<N  } 
and 

T  -  { (i,i,i-l,i):0<isN  j  . 

4.  THE  DECOMPOSITION  METHOD 

The  marked  job  method  of  Section  3  is  applicable  to  passage  times  in  the  general  sense,  i.e.. 
whether  or  not  the  passage  time  is  a  complete  circuit  in  the  network.  For  general  characteristics  of 
passage  times  which  are  complete  circuits,  no  other  method  for  obtaining  confidence  intervals  from  a 
single  simulation  run  is  known  to  the  authors.  Note  that  since  only  the  observed  passage  times  for  the 
marked  job  enter  into  the  construction  of  point  and  interval  estimates,  there  may  be  some  loss  of 
statistical  efficiency  as  the  price  for  obtaining  confidence  intervals. 


II 


In  this  section  we  concentrate  on  passage  times  through  a  subnetwork  of  a  given  closed  network 
of  queues,  i.e.,  passage  times  which  are  not  complete  circuits  in  the  network.  For  this  class  of  passage 
times  we  develop  an  estimation  method  in  which  observed  passage  times  for  all  the  jobs  enter  into  the 
construction  of  point  estimates  and  confidence  intervals.  We  consider  closed  networks  of  queues  and 
passage  times  as  in  Section  2.  but  we  make  the  additional  assumption  with  respect  to  the  sets  S  and  T  of 
(5)  that  S~T-o;  this  effectively  rules  out  passage  times  which  are  complete  circuits. 

The  basis  of  the  decomposition  method  for  estimation  of  passage  times  through  a  subnetwork  is 
simulation  of  the  network  in  random  blocks  defined  by  the  terminations  of  certain  passage  times.  The 
distinguished  passage  times  arc  those  that  (i)  terminate  when  no  other  passage  times  are  underway,  and 
(ii)  leave  a  fixed  configuration  of  the  job  stack  defined  by  (1).  These  terminations  serve  to  decompose 
the  sequence  of  passage  times  for  all  of  the  jobs  into  independent  and  identically  distributed  blocks. 

We  denote  by  { Pn°:n2l  j  the  sequence  of  passage  times  (irrespective  of  job  identity)  enumerated 
in  order  of  passage  time  start.  As  before,  we  let  f  be  a  real-valued  function  with  domain  R_,  and  the 
goal  of  the  simulation  is  the  estimation  of 

r°(f)  -  E  j  f(P°)  J  ,  (9) 

where  Pn°-*P°.  Note  that  P°-P.  the  limiting  passage  time  for  (any)  marked  job. 

Algorithm  2.  Decomposition  method 

1.  Select  a  configuration  z°  of  the  job  stack  at  which  a  passage  time  terminates  and  there  are  no  other 
passage  times  underway  .  Begin  the  simulation  with  this  configuration  of  the  job  stack  at  time  0. 

2.  Carry  out  the  simulation  for  a  fixed  number  n  of  blocks  defined  by  the  successive  terminations  of 
passage  times  irrespective  of  job  identity  which  leave  the  job  stack  in  the  (fixed)  configuration  z°. 


3.  In  each  block,  measure  the  passage  times  for  all  the  jobs  and  record  these  along  with  the  number  of 
passage  times  observed  in  the  block. 
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4.  Denote  the  number  of  passage  times  observed  in  the  mth  block  b>  Km°  and  compute  Ym°(f).  the  sum 
of  the  quantities  f(Pj°)  for  the  passage  times  Pj0  in  the  mth  block. 

5.  Based  on  n  blocks,  form  point  and  interval  estimates  for  r°(f)  from  the  quantities 
{ i  Ym°(f).Km°):l£msn  J  as  in  the  standard  regenerative  method. 


The  underlying  stochastic  structure 

We  begin  by  labelling  the  jobs  from  1  to  N.  and  for  i— 1,2 . N,  denote  by  N‘U)  the  position  of 

job  i  in  the  linear  job  stack  at  time  t.  Then  in  terms  of  the  vector  Z(t)  of  ( 1 ),  for  t^O  we  set 

X°(t)  -  (Z(t).N'(t) . NN(t)).  (10) 

The  process  X°- { X°(t):t20 )  is  an  irreducible,  postive  recurrent  continuous  time  Markov  chain  with 
finite  state  space  E°.  Next  we  let  L°(t)  denote  the  last  state  visited  by  by  the  Markov  chain  X°  before 
jumping  to  X°(t).  and  for  tsO  define 


V»(t)  -  (L°(t),X°(t))  . 


(ID 


The  process  V°-  j  V°(t):tss0  J  is  the  fundamental  stochastic  process  of  the  simulation.  Since  X°  is  an 
irreducible,  positive  recurrent  continuous  time  Markov  chain,  so  is  V°.  We  denote  the  state  space  of  V° 
by  F°  and  define  two  subsets  S°  and  T°  of  according  to 

S°  -  { (z.n  | nvz'.n|' . nN')«F°:  for  some  k.  (z.nk)t  A,  and  (z'.nk')«  A;j 


and 


T°  -  {(z.n, . nN,z'.n,' . nN')«F°:  for  some  k,  (z.nk)«  B,  and  (z.nk')«  B;  J  .  (12) 

The  entrances  of  V°  to  the  set  S°  correspond  to  the  starts  of  passage  times  (irrespective  of  job  identity) 
and  the  entrances  of  V°  to  the  set  T°  correspond  to  the  terminations.  Thus  from  a  simulation  of  the 
process  V°.  it  is  possible  to  measure  the  passage  times  for  all  the  jobs. 


Now  consider  [Pn°:nal  |  the  sequence  of  passage  times  (irrespective  of  job  identity),  enumerat¬ 
ed  in  order  of  passage  time  start.  Recall  that  the  goal  of  the  simulation  is  estimation  of 

r°(f)  -  E  { f(P°) }  . 

where  Pn°~P°.  and  f  is  a  real-valued  measurable  function  with  domain  R.. 
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^  e  cam  out  the  simulation  in  random  blocks  of  the  process  V®  defined  by  the  successive 
entrances  of  to  a  fixed  set  of  states  U°.  Entrances  of  V°  to  the  set  U°  correspond  to  the  terminations 
of  passage  times  (irrespective  of  job  identity)  which  occur  when  no  other  passage  times  are  underway, 
and  which  leave  a  fixed  configuration  (z°)  of  the  job  stack:  see  [10]  for  a  formal  definition  of  the  set  L'°, 
For  convenience,  we  assume  that  V°(0)«L0.  Denoting  by  j*ym°:m2il  [  the  lengths  in  discrete  time  units 
of  the  successive  blocks  (returns  to  the  set  U°)  for  [  Vn°:niOJ  the  embedded  jump  chain  associated  with 
'°-  we  define  o0°-0  and  <5m°— ri°-...-,-ym0.  matl.  (Note  that  the  successive  entrances  to  the  set  U°  are 
not  regeneration  points  for  the  process  V°).  The  number  of  passage  times  K,°  in  the  first  block  of  the 
process  V°  is 


and  we  denote  the  analogous  quantity  in  the  mth  block  of  V°  by  Km°.  Note  that  within  each  block  of  V° 
defined  by  entrances  to  the  set  U°.  at  least  one  passage  time  starts  and  terminates.  Next,  we  let  Y  °(f) 
be  the  sum  of  the  quantities  f(Pj°)  over  the  passage  times  Pj°  in  the  mth  block  of  V°.  e.g.. 

K,° 

Ym°(f)  -  I  f(P»)  . 
j-i 

The  key  observation  is  that  the  sequence  of  pairs  of  random  variables 

[ (Ym°(0.Km°):maI }  (13) 

are  i.i.d.  For  the  function  f  appearing  in  (9).  let  D(f)  denote  the  set  of  discontinuities  of  f.  Assuming 
that  P )  P°  t  D(  0  }  “0.  it  follows  that  as  n— co, 

f(Pn°)—f(P°)  . 

Finally,  provided  that  P{P«D(f)]-0  and  E(  |  f(P°)|  }<oc,  standard  arguments  (cf.  [8])  establish  the 
ratio  formula 


r°(f)  -  E{  Y ,°(f)  ]  /E{  K,°)  . 
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Given  this  ratio  formula  and  the  fact  that  the  random  variables  in  (13)  are  i.i.d..  the  regenerative 
method  provides  from  a  fixed  number  of  blocks  a  strongly  consistent  point  estimate  and  an  associated 
confidence  interval. 

5.  EFFICIENCY  OF  SIMULATION 

For  mean  passage  times,  theoretical  values  can  be  obtained  [10]  for  variance  constants  entering 
into  the  central  limit  theorems  used  in  previous  sections  to  obtain  confidence  intervals  for  passage  time 
characteristics.  These  calculations  provide  a  firm  basis  for  an  assessment  of  the  statistical  efficiency  of 
the  two  methods. 

For  a  comparison  of  the  statistical  efficiency  of  the  marked  job  and  decomposition  methods 
where  both  apply,  it  is  convenient  to  have  a  central  limit  theorem  comparable  to  (81  but  in  terms  of 
simulation  time,  t.  rather  than  number  of  cycles,  n.  Let  m(t)  be  the  number  of  passage  times  completed 
by  time  t.  i.e..  in  the  interval  (0,t].  Then  as  t-oo, 

t1 ' 2[ {  m(t)  }'*(  "I* *  f(Pj))-r(0]/[(E {a, } )*''2<r(0/E {  M, } ]  -  N(0,1)  . 
i-l 

where  Eja,  ]  is  the  expected  length  of  a  0-cycle  in  the  continuous  time  process  V.  It  follows  that  the 
quantity 

e  -  (E[a,  j 

is  the  appropriate  measure  of  statistical  efficiency  for  the  marked  job  method  and  is  independent  of  the 
state  0«S  selected  to  form  cycles.  (For  mean  passage  times,  the  function  f  is  the  identity  function  and 
<r-  denotes  the  corresponding  variance  constant.)  Similarly,  the  quantity 

eOMElVn'/WE^O}  . 

where  a,0  is  the  length  of  a  block  in  V°.  is  the  appropriate  measure  of  statistical  efficiency  for  the 
decomposition  method. 


Table  I  gives  theoretical  values  for  simulation  by  the  decomposition  method  for  r°-E  [  P°  |  in  the 
network  of  queues  of  Figure  2.  Results  are  given  for  N-i  to  6  jobs;  A,  and  A-  are  the  rate  parameters 
for  the  exponential  service  times  at  centers  1  and  2,  respectively,  and  p  is  the  probability  of  feedback  to 
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center  I.  We  hold  the  values  of  X,-1.0  and  p-0.75  fixed,  but  vary  X,.  For  N-2  jobs  (with  p-0.75. 
X,-1.0  and  X:-0.5).  the  value  of  e°  which  measures  the  statistical  efficiency  of  the  decomposition 
method  is  16.546.  The  corresponding  value  for  the  marked  job  method  is  20.890.  Thus,  for  these 
parameter  values  the  decomposition  method  is  approximately  21  percent  more  efficient  than  the  marked 
job  method.  For  N-4  jobs,  the  decomposition  method  is  41  percent  more  efficient.  Table  2  gives  a 
comparison  of  the  relative  efficiency  (e/e°)  of  the  marked  joD  and  decomposition  methods  for  the  same 
sets  of  parameter  values  as  in  Table  1. 

6  CONCLUDING  REMARKS 

W'e  have  limited  the  discussion  of  estimation  method  for  passage  times  to  networks  of  queues 
having  service  centers  which  operate  as  a  single  server.  The  generalization  to  networks  having  multiple 
server  service  centers,  however,  is  straightforward.  To  estimate  passage  times  it  is  sufficient  to 
incorporate  into  the  linear  job  stack  (and  state  vector  definition)  information  carrying  the  class  of  job 
being  serviced  by  each  of  the  servers  at  a  mutiple  server  service  center. 

Estimation  methods  for  passage  times  in  finite  capacity,  open  networks  of  queues  are  available. 
In  [8],  two  formulations  of  the  finite  capacity  constraint  are  considered  and  particular  stochastic  point 
processes  associated  with  a  Markov  renewal  process  generate  arrivals  to  the  networks.  These  so-called 
Markov  arrival  processes  facilitate  the  incorporation  into  computer  and  communication  system  models  of 
departures  from  a  Poisson  arrival  pattern.  The  basis  for  the  estimation  of  passage  times  is  the  measure¬ 
ment  of  passage  times  for  a  sequence  of  marked  jobs;  these  are  typical  jobs  in  the  sense  that  the 
sequence  of  passage  times  for  the  marked  jobs  coverges  in  distribution  to  the  same  random  variable  as 
do  the  passage  times  for  all  jobs. 

The  marked  job  method  of  Section  3  applies  equally  well  to  networks  with  stochastically 
nonidentical  jobs.  In  such  networks,  jobs  of  each  type  have  their  own  routing  structure  and  service 
requirements,  and  in  the  case  of  finite  capacity,  open  networks,  (independent)  arrival  processes. 
Methods  for  the  estimation  of  joint  characteristics  of  passage  times  over  the  several  job  types  are 
discussed  in  [9]. 
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TABLE  1 

Statistical  efficiency  of  the  decomposition  method 


P" 

V 

^2' 

0.75 

-1.0 

-0.125 

p-0.75 

X,-l  0 
X-,-0.25 

p-0.75 

X,  —  1 .0 

X-,-0.5 

N 

E  { P) 

e° 

E I P  | 

e° 

E  { P I 

eO 

1 

4.000 

13.856 

4.000 

1 1.314 

4.000 

9.798 

2 

5.333 

19.956 

6.000 

17.664 

6.667 

16.546 

3 

6.286 

27.380 

8.000 

26.128 

9.714 

25.035 

4 

6.933 

35.189 

10.000 

36.606 

13.067 

34.491 

5 

7.355 

42.597 

12.000 

49.107 

16.645 

44.296 

6 

7.619 

49.068 

14.000 

63.645 

20.381 

54.021 
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TABLE  2 

Relative  efficiency  of  marked  job  and  decomposition  methods 

p-0.75 

p-0.75 

p-0.75 

A,-1.0 

A, -1.0 

A ,  —  1 .0 

A:-0.I25 

A,-0.25 

A:-0.5 

N 

1 

1.000 

1.000 

1.000 

2 

1.190 

1.189 

1.211 

3 

1.207 

1.224 

1.319 

4 

1.190 

1.209 

1.408 

5 

1.176 

1.186 

1.499 

6 

1.394 

1.162 

1.597 
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Figure  1.  Exponential  stage  reoresentation 


Figure  2.  Closed  network  of  queues 
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